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LINEAR INTERPOLATION OF FOUR-DIMENSIONAL 
TABULATED DATA FOR COMPUTERS WITH SINGLE 
SUBSCRIPTED VARIABLE CAPABILITY 

INTRODUCTION 

At present several numerical methods exist for interpolation of 
tabulated data of functions of several variables, notable of which are the 
Gregory-Newton, Gauss, and Lagrangian. These methods depend, however, 
upon the calculation of finite differences or coefficients to fitted polynomials. 

If the tabulated data are approximately linear, the above methods are very 
accurate but the computation time involved is excessive, especially for digital 
computer applications. Linear interpolation has approximately the same 
accuracy but with much less computation time. 

With the purpose of faster computation, a Fortran digital computer sub- 
program was developed to linearly interpolate tabulated data of four or fewer 
dimensions. As a test of the program' s proficiency and accuracy, it was 
used to interpolate approximately linear aerodynamic data which were tabulated 
as a function of four variables. The results of the linear interpolations were 

compared with those of Lagrangian interpolations of the same data. The 
answers varied only in the fourth decimal place; however, the linear routine 
extracted approximately 1000 values in the time it took the Lagrangian routine 
to extract 400 values. The slight inaccuracy of the linear method was offset 
by its inherent speed. 


MATHEMATICAL MODEL 

The first step in establishing a mathematical model for the linear 
interpolation of a function Q of the four independent variables X, Y, Z, and W 
in tabulated form would be to express Q as a subscripted variable of the 
fourth dimension: 

Q = Q(X(I), Y(J), Z(K), W(L)) (1) 

with a single subscripted variable for each of the independent variables, 
where I, J, K, and L are the number of elements in each array, respectively. 



Linear interpolation could then be accomplished easily by finding the location of 
the desired point Q (XS, YS, ZS, WS) from the independent variable arrays. 

This method is fine for large computers, but it fails for many smaller machines 
which allow a maximum of only three subscripts. To permit utilization of the 
method on any computer with one -dimensional subscripted variable capability, 
the four-dimensional Q array may be mapped into a one -dimensional array 
Q* (N) (the * notation will be dropped for brevity in the following discussion) . 

The mapping is as follows, selecting a (3 x 3 x 3 x 3) array as an 
example: 

Q (1) - Q (1, 1, i, i) 

Q (2) = Q ( 1 , 1, 1, 2) 

Q (3) = Q (1 , 1, 1, 3) 

Q (4) = Q(l, 1, 2, i) 

Q (5) = Q(l, i, 2, 2) 

Q (6) = Q(l, 1, 2, 3) 

Q (7) = Q (1 , i, 3, 1) 

(2) 

Q (8) - Q (1 , 1, 3, 2) 

Q (9) - Q(i, 1, 3, 3) 

Q (10)= Q (1, 2, 1, 1) 

Q ( 1 1) - Q (1 , 2, 1, 2) 


Q (N')= Q(I\ J\ K\ L’> 

N' = (!’ - 1) * J • K • L + (J' - 1) • K * L + (K' - 1) • L + L' . (3) 
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For the (3 x 3 x 3 x 3) example there will be I • J • K • L or 81 elements in the 
Q(N) array. Figure 1 depicts the complete mapping for the (3 x 3 x 3 x 3) 
example. However, because of the form of certain mapping equations [see 
equations (5) ], the number of elements in Q must be increased by 
L • (K + i) + 1 . These added terms are for working storage, and their values 
must be initially zero. Therefore, the number of elements in Q must be 


N = 

i • 

J * 

K 

• L + 

L 

(K + i) 

+ 1 . 
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Figure 1. Four-dimensional (3 x 3 x 3 x 3) array 
mapped into a one-dimensional array. 
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Associated with the one-dimensional Q array will still be a one- 
dimensional array for each of the independent variables X, Y, Z, andW. These 
arrays contain the data points along the respective independent variable ranges, 
e.g. , X may vary from 10 to 100 in increments of 10, thus giving 10 values for 
the X array. The elements of these arrays must be in nondecreasing order. 

METHOD OF LINEAR INTERPOLATION FOR 
FOUR DIMENSIONS 

After the dependent and independent variable arrays have been estab- 
lished, the dependent variable Q may be determined for a given point with inde- 
pendent variables (XS, YS, ZS, VVS) by linear interpolation. 

Initially, the location of the point's independent variable coordinates 
within the respective independent variable arrays must be ascertained. For 
example, XS obviously lies between some X(n) and X(II + 1); also YS lies bet- 
ween some Y (JJ) and Y (JJ + 1) , etc. The indices of prime interest are II, JJ, 
KK, and LL. If any of the point's coordinates happen to coincide with the last 
entry in the respective array, that corresponding index is equated to the array 
dimension, e.g., 11 = I. Also, if any of the point's coordinates lie outside the 
range of their respective arrays, these variables are assumed to be the end- 
points . 

Again, referring to the (3 x 3 x 3 x 3) array in Figure 1, the linear 
interpolation method may be illustrated with, sample values for (II, JJ, KK, 

LL) of (1, 2, 1, 2 ) . These known indices facilitate the determination of the 
appropriate values of Q, which must be interpolated in each of the four dimen- 
sions. As noted in Figure 1, these particular elements of the Q array are 
subscripted I x , I 2 , I 3 , I 4 , I’p I’ 2 , I' 3 , I’ 4 , and are obtained from the follow- 
ing mapping equations: 


I t = (II - 1 ) • J * K * L + (JJ - i) • K • L + (KK - 1) • L + LL 


1 2 = II + L 

1 3 = I t + K * L 


(5) 


I 4 = I 3 + L 
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The primed subscripts are obtained by adding J • K • L to the preceding unprim 
ed subscripts, respectively. The equation for I 4 forces the Q storage require- 
ment increase referred to above. 

A schematic of the necessary linear interpolations in each of the four 
dimensions is shown in Figure 2. Interpolation of the aforementioned sub- 
scripted elements of Q in the W dimension yields the values T 1# T 2 , S l5 T'j, 

T' 2 , S'j, S' 2 . These values are then interpolated in the Z dimension to pro- 
duce Uj, U 2 , U'j, U' 2 , which are interpolated to obtain Vj and V 2 in the Y 
dimension. Q(XS, YS, ZS, WS) is then obtained from the X dimension inter- 
polation of V t and V 2 . The linear interpolations are of the form 


Q = Vj + 


XS - X(I I) 

X (I I + 1 ) - X(I I) 


( 6 ) 


The required value of Q is thus easily obtained once the appropriate Q sub- 
scripts are known. 


METHOD OF LINEAR INTERPOLATION IN 
FEWER THAN FOUR DIMENSIONS 


The previous equations pertaining to linear interpolation of four- 
dimensional tabulated data may also be used with 1, 2, or 3 dimensions. The 
only necessary requirements for three dimensions are that L and LL be equated 
to one, W(l) be equated to zero, and W(2) be equated to a positive real number 

In two dimensions the three-dimensional requirements must be met in 
addition to similar requirements on the Z dimension parameters: K = 1, 

KK = 1, Z(i) = 0, and Z(2) > 0, with additional simplifications: 


Is =Ii 


I 4 = I 2 


(7) 


For one dimension the two-dimensional requirements must be met plus 
similar requirements on the Y parameters: J = 1, JJ = 1, Y(l) = 0, Y (2) > 0. 
There is no need now for the primed indices on Q since the interpolation 
equations need hold for only the unprimed indices. 
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COMPUTER SUBPROGRAM 


Two Fortran computer subprograms were written for linear interpola- 
tion of tabulated data of four or fewer dimensions. Listings of these routines 
are given in Appendix A and B . The routine referred to in Appendix A requires 
no more than single-subscript capability, whereas that in Appendix B requires 
double -subscript capability. The computer word storage allocations for the 
respective routines are Routine A, 727, and Routine B, 718. Routine B is pre- 
ferred because of its smaller storage allocation and slightly faster computation 
time. However, it is somewhat more complicated in the expression of its argu- 
ments . 


The included comment cards at the beginning of each routine should be 
sufficient to enable easy use of the methods based on the mathematical model 
presented herein. 
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APPENDIX A 


FOUR-DIMENSIONAL LINEAR INTERPOLATION 
SUBPROGRAM FOR COMPUTERS WITH SINGLE 
SUBSCRIPTED VARIABLE CAPABILITY 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FUNCTION TLU4D ( N, J, K , L , X , V , Z , U , 0, XS , Y S , Z S , WS ) A 1 

dimension x< 1), Y(l), Z 1 1 ) , w<l), 0(1), V < 2 ) A 2 

A 3 
A A 

FOUR-DIMENSIONAL TaRLE LOOK-UP routine A 5 

A 6 

0 IS A FUNCTION OF THE FOUR VARIABLES X,Y,Z. AND W, IN THE A 7 

table, a one dimensional array is entered for each variable and a s 

0. the dimensions of the arrays are specified as...' a 9 

A 10 

N number of elements IN X ARRAY A 11 

J number of elements in y ARRAY A 12 

K NUMBER OF ELEMENTS in l 'array A 13 

L NUMBER of ELEMENTS IN W ARRAY A 14 

A 15 

the VALUES OF a IN THE TABLE ENTRY array ARE EXPRESSED AS in A 16 

THIS Example of a 3. 3. 3, 3 ARRAY,,. ' A 17 

A 18 

0(1) S QCl.1.1,1) A . 19 

Q ( 2 ) = 0(1. 1.1, 2) A 20 

0(3) = 0(1. 1.1.3) A 21 

0(4) = 0(1. 1.2,1) A 22 

Q ( 5 ) = 0(1. 1.2. 2) A 23 

0(6) * 011,1,2,3) A 24 

0(7) *0(1. 1.3.1) A 25 

0(d) =0(1, 1.3. 2) A 26 

0(9) = 00,1.3,3) A 27 

0(10) = 0(1, 2.1,1) AND SO LN, A 28 

A 29 

THE DIMENSION OF 0 MUST BE N*J»K»L * L*(K+1)+1 WHERE THE A 30 

L* ( K*1 > *1 TERM IS FOR WORKING STORAGE, THE ELEMENTS OE THE A 31 

WORKING STORAGE AREA OF TH£ ARRAY MUST PE INITIALLY ZEROED OUT, A 32 

A 33 

given values for the four variables (xs,ys,zs,ws> THt routine a 34 

LINEARLY INTERPOLATES for the value of OIXS. Ys, ZS. WS j , THE A 35 

IF any OF the FOUR VARIABLES XS.YS.ZS.WS LIE OUTSIDE the A 36 

range of their respective arrays, these variables ARt assumed a 37 

TO BE THE ENDPOINTS, * A 38 

ANSWER BEING EXPRF-SSED as TLU4D. A 39 


A 40 

THE ROUTINE is WRITTEN FOR FOUR DIMENSIONAL USAGE, BUT MAY M t A 41 

USED for 1,2, and 3 DIMENSIONS. FOR A THREE-DIMENSIONAL TABLE A 42 

L MUST BE SET EQUAL TO 1 AND A DUMMY VARIABLE D(K> MLST BE SET A 43 

FOR THE W ARRAY AS 0 ( 1 ) = 0 , 0 AND D(2)= ANY PCSI T I VE' VALUE A 44 

greater than o. for a two-dimensional table k and l must be a 45 

EOUAL TO 1 WHILE The dummy ARRAY replaces BOTH the w and Z A 46 

arrays, foh a one-dimensional Table u,k,l must all be set a 47 

EQUAL TO 1 AND THE DUMMY ARRARY REPLACES the Y ,1, ANl W ARRAYS. A 48 

A 49 
A 50 
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1=1 A 51 

1U IF (X( I ) - X S > 20,20.50 A 52 

20 IF < !-M) 30,40,40 A 53 

30 1=1+1 A 54 

00 TO 10 A 55 

4 1) I I = N A 55 

no TO 60 A 57 

50 I 1 = 1-1 A 58 

60 IF (J-l) 70,70,80 A 59 

7 0 JJ = 1 A 60 

KK = 1 A 61 

LL=1 A 62 

|<C' TO 30 0 A 63 

C A 64 

an i=i a 65 

VO IF (Y(I)-YS) 100,100,130 A 66 

100 IF (I-J) 110,120,120 A 67 

110 1=1+1 A 68 

GO TO 90 A 69 

l2l) JJ = j A 70 

no TO 140 A 71 

1 3 1) JJ»I-1 A 72 

140 IF (K-l) 150,150,160 A 73 

l5n K K = 1 A 74 

LL=1 A 75 

HO TO 300 A 76 

C A 77 

16(1 1=1 A 78 

1711 IF ( Z ( I ) ■ 2 S ) 180.180,210 A 79 

16 ll IF ( I - K ) 190,200,200 A 80 

1 V0 1=1+1 A 81 

GO TO 170 A 62 

200 KK =K A 83 

GQ TO 220 A 84 

210 KK=I-1 A 65 

220 IF ( L - 1 ) 230.23C.P40 A 66 

230 tL=i A 87 

G(j Tn 300 A 88 

C A 89 

240 '1 = 1 A 90 

2 5 11 IF (W(I)-WS) 26n,?6o.29c A 9l 

260 IF (I-l) 270,280,260 A 92 

270 1=1*1 A 93 

no TO 250 A 94 

260 L L * L A 95 

GlJ TO 300 A 96 

2V0 LL= I *1 ' A 97 

C A 98 

300 MN=J+K+L A 99 

IF (II) 320,310,320 A 100 

310 11=1 A 101 

320 IF <JJ> 340,330,340 A 102 

330 JJ=1 A 103 

340 IF (KK) 360,350,360 A 104 

350 KK = 1 A 105 
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360 

IF <t_D 380,370,380 

A 

106 

370 

LL = 1 

A 

107 

360 

Il=(!l-1>*NN*( JJ-1>*K*L*(KK-1)*L*LL 

A 

108 


1 2= 1 1*L 

A 

109 


IF (K-l) 390,390,400 

A 

110 

390 

13=11 

A 

111 


14=12 

A 

112 


<50 TO 410 

A 

113 

400 

I3=U*K*L 

A 

114 


14=13*1 

A 

115 

C 


A 

116 

410 

00 440 1=1, 2,1 

A 

117 


WlB(WS-WtLt) )/(W<LL*l)-W(LU) ) 

A 

118 


Tl = 0( I1)*W1*(Q( I1*1)«Q( 11) ) 

A 

119 


T2«0( 12 )*M1* (0(12*1 )»0( 12)) 

A 

120 


S1 = Q( I3)*W1*(q( I 3*l)-Q( 13) > 

A 

121 


S2=0( !4,*wl*(Q< I4*1)*Q(I4)) 

A 

122 


II1 = T1*(2S«Z(KK) )«(T2»T1)/(Z<KK*1 > - 2 < KK > ) 

A 

123 


U2=S2*(ZS*Z(KK) >*(S2*S1>'<Z<KK*1)-Z<KK) ) 

A 

124 


V< I >=U1*<VS-Y<JJ) )*(U2-U1)/(Y(JJ*1)*Y(JJ) ) 

A 

125 


IF (J-l) 420,420,430 

A 

126 

420 

V < 2 ) =u2 

A 

127 


GO TO 450 

A 

128 

430 

11= 1 1*NN 

A 

129 


12= 12 + NIN 

A 

130 


13=13* NN 

A 

131 


14=14* NN 

A 

132 

440 

CONTINUE 

A 

133 

450 

TLU40 = V(1 )*(XS*X( I I ) )*(V(2)”V<1) )/(X( I I *1 ) -X ( l I ) ) 

A 

134 


RETURN 

A 

135 


END 

A 

136 
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APPENDIX B 


FOUR-DIMENSIONAL LINEAR INTERPOLATION 
SUBPROGRAM FOR COMPUTERS WITH DOUBLE 
SUBSCRIPTED VARIABLE CAPABILITY 



oo nnoooooonooooononnno noooorioocjoonooooooooooonno 


FUNCTION TLU4D <M,X,Q,XS> B 1 

DIMENSION 0(1), V ( 2 ) , P( 5), XS(1), M ( 1 ) , L<4), X(l,l) B 2 

B 3 
B 4 

four-dimensional table look-up routine b 5 

B 6 
B 7 

q is a function of the four variables x.y.z.w and h is a one- b e 

DIMENSIONAL array with four elements WHICH ARE THE NUMBER OF B 9 

ELEMENTS in the independent variable ARRAYS ' B 10 

B 11 

N s Mil) = NUMBER OF ELEMENTS IN X ARRAY B 12 

J = M(2> = NUMBER OF ELEMENTS -IN Y ARRAY B 13 

K = H ( 3 ) * NUMBER OF ELEMENTS IN Z ARRAY B 14 

L = M<4) = NUMBER OF ELEMENTS IN W ARRAY B 15 

B 16 

WHERE ... B 17 

X ( 1 , N ) CORRESPONDS TO THE X ARRAY 8 18 

X(2,N) CORRESPONDS TO THE Y ARRAY B 19 

X ( 3, N ) CORRESPONDS TO THE Z ARRAY B 20 

X ( 4 , N I CORRESPONDS TO THE W ARRAY B 21 

B 22 

AND ALSO ... 8 23 

XS(1) CORRESPONDS TO THE POINT COORDINATE XS B 24 

XS ( 2 ) CORRESPONDS TO THE POINT COORDINATE YS B 25 

X S < 3 ) CORRESPONDS TO THE POINT COORDINATE ZS B 26 

XS<4) CORRESPONDS TO THE POINT COORDINATE MS 9 27 

B 28 
B 29 

THE DIMENSION OF 0 MUST BE N*J*K*L + L*(«*1>*1 WHERE THE B 30 

L* ( K+l ) +i TERM IS FOR WORKING STORAGE. THE FLEMENTS OF THE B 31 

WORKING STORAGE AREA OF THE ARRAY MUST RE INITIALLY ZEROED OUT, B 32 

B 33 

given values For the four variables (XS. ys.zs.wsj the routine b 34 

linearly interpolates for the value of o<xs,ys,zs,ws>. the b 35 

answer being EXPRESSED AS TLU4D, B 36 

IF ANY OF T H£ FOuR VARIABLES XS.YS.ZS.wS LIE CutSIDE the B 37 

range of their respective arrays, these variables are assumed b 38 

TO BE the endpoints, B 39 

B 40 

THE ROUTINE IS WRITTEN FOR four dimensional usage, BUT MAY BE B 41 

USED for 1.2 , and 3 DIMENSIONS, FOR A THREE-DIMENSIONAL TABLE B 42 

L must BE SET EQUAL TO 1 AND A DUMMY VARIABLE D ( K ) MUST BE SET B 43 

FOr THE W ARRAY AS D<1>=0,0 AND C(2)= ANY POSITIVE VALUE B 44 

GREATER THaN 9, FOR A TWO-DIMENSIONAL TABLE K AND L MUST BE 8 45 

EQUAL TO 1 WHILE The DUMMY ARRAY REPLACES both THE w and Z B 46 

ARRAYS. FOR A ONE » D I MENS I ONA L TABLE J,K«L MUST ALL BE SET B 47 

EQUAL TO 1 AND THE DUMMY ARRARY REPLACES THE Y.Z, ANC W ARRAYS. B 48 

B 49 
B 50 
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DO 100 J=l,4,l 

B 

51 


I=M< J) 

B 

52 


00 in K = 1 , [ , 1 

B 

53 

in 

P(K)=X(J,K) 

B 

54 


PS=XS ( J ) 

B 

55 


IF (M(J)-l) 80,80,20 

B 

56 

20 

1=1 

B 

57 

3 0 

IF (P(I)-PS) 40,40,70 

B 

58 

40 

IF ( I-M(J) ) 50,60,60 

B 

59 

5d 

1 = 1*1 

B 

60 


GO TO 3q 

B 

61 

60 

L ( J) =M( J) 

8 

62 


GO TO 100 

B 

63 

70 

U J)=I-1 

B 

64 


GO TO 100 

B 

65 

60 

00 90 K = J , 4 , 1 

B 

66 

90 

1- ( K ) = 1 

B 

67 


GO TO 130 

B 

68 

100 

CONTINUE 

B 

69 


DO 120 K=1 , 4 , 1 

B 

70 


IF < L < « ) ) 120,110,120 

B 

71 

11(1 

L ( K ) = 1 

8 

72 

120 

CONTINUE 

B 

73 

1 3 n 

N = M(2)*M(3)*iM(4) 

B 

74 


Il=a(l)*l>*fJ*(L(2)«l)*M(3)*M(4)*(L(3>-l)*M(4)*L(4) 

B 

75 


1 2 = 1 1 + M ( 4 ) 

B 

76 


IF <M(3)-D 140,140,150 

B 

77 

140 

13=11 

B 

78 


14=12 

B 

79 


GO TO 160 

B 

80 

1 5 0 

I 3= 1 1*M( 3 ) *H( 4 ) 

B 

81 


I 4 = I 3 ♦ M < 4 ) 

B 

82 

160 

DO 190 K=1 , 2 , 1 

B 

83 


I I =U ( 1 ) 

B 

84 


JJ=L<2> 

B 

85 


*K=L(3) 

B 

86 


LL=U<4) 

B 

87 


Wl=(XS<4)*X(4,tL))/(X(4,LL*l)-X(4,LL)) 

B 

88 


Tl = 0(U)*W.l*(Q(Il*l)-0<Il)) 

B 

89 


T2 = Q( I 2 ) *W1* ( o ( I 2*1 ) =Q ( 12) ) 

B 

90 


S1 = 0( 13)+W1*(0( 13*1)»0( 13) ) 

B 

91 


S2 s Ql I 4 ) *wi* ( Q ( I4 *i)sQ(J4) ) 

B 

92 


Ul=Tl*IXS(3)-x(3,KK) )*<T2-T1)/(X(3,KK*1)*X(3,KK) ) 

B 

93 


U2=S2+(XS(3)-X(3,KK) )*(S2-S1)/<X(3,KK*1)-X(3,KK)) 

B 

94 


V(K)=ul*(XS(2)-X(2,JJ) )*(U2-U1)/(X(2, JJ*1 >-X.( 2 > JJ) ) 

B 

95 


IF <M(2)-1) 170,170,180 

B 

96 

1 7 o 

V ( 2 ) =U2 

B 

97 


GO TO 200 

B 

98 

160 

U = U*N 

B 

99 


Z2. 

♦ 

cv 

11 

CM 

B 

100 


I 3 = I 3 ♦ N 

H 

101 


I 4= 14*N 

B 

102 

19H 

CONTINUE 

B 

103 

200 

TLU4D=V(1)*(XS(1)-X(1,II))*(V(P)-V(1))/(X(1,I1*1)-X(1,U)> 

B 

104 


RETURN 

B 

105 


END 

B 

106 
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